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Abstract 

The goal of this tutorial is to promote interest in the study of 
random Boolean networks (RBNs). These can be very inter- 
esting models, since one does not have to assume any func- 
tionality or particular connectivity of the networks to study 
their generic properties. Like this, RBNs have been used 
for exploring the configurations where life could emerge. 
The fact that RBNs are a generalization of cellular automata 
makes their research a very important topic. 

The tutorial, intended for a broad audience, presents the state 
of the art in RBNs, spanning over several lines of research 
carried out by different groups. We focus on research done 
within artificial life, as we cannot exhaust the abundant re- 
search done over the decades related to RBNs. 

Introduction 

Random Boolean networks (RBNs) were originally devel- 
oped by Stuart Kauffman as a model of genetic regulatory 
networks ( Kauffman, 1969 Kauffman, 1993 1. They are also 
known as N — K models, or Kauffman networks. RBNs are 
generic, because one does not assume any particular func- 
tionality or connectivity of the nodes: these are generated 
randomly. This is a useful approach if the specific structure 
and/or function of a system are very complex and unknown. 
The generic properties found in the model can be then ap- 
plied to the particular system, in order to attempt to unveil 
its mechanisms. 

Mathematical and computational modelling of genetic 
regulatory networks promises to uncover the fundamental 
principles of living systems in an integrative and holistic 
manner. It also paves the way toward the development of 
systematic approaches for effective therapeutic intervention 
in disease ( Sh mulevich et al., 2002| . Single-gene studies are 
very limited for such intertwined networks. 

Even when there is a Boolean simplification, many sys- 
tems can be studied with near-binary states. This is because 
the behaviour of many systems is determined by thresholds, 
such as the ones determined by firing potentials of synapses 
in neurons, or activation potentials of chemical reactions in 
metabolic networks. 



Furthermore, random Boolean networks have been ap- 
plied and used as models in many different areas, such 
as evolutionary theory, mathematics, sociology, neural net- 
works, robotics, and music generation. 

In the next section, we will review the classic RBN model, 
its three characteristic phases (ordered, chaotic, and criti- 
cal, some explorations of the model, alternatives, and ex- 
tensions. Next we will review the effect of the updating 
scheme in RBNs: syncrhonous-asyncrhonous, determiistic- 
non-deterministic. We mention briefly some applications of 
RBNs, tools availabe for their study, and future lines of re- 
search. 

Classical Model 

Kauffman proposed the original RBN model, supporting the 
hypothesis that living organisms could be constructed from 
random elements, without the need of precisely programmed 
elements (Kauffman, 1969). Certain types of RBNs are very 
robust, and have many analogies to living organisms. 

A RBN consists of N nodes, which can take values of zero 
or one (Boolean). The state (zero or one) of each node is de- 
termined by K connections coming from other (or the same) 
nodes. The connections are wired randomly, but remain 
fixed during the dynamics of the network, i.e. "quenched". 
The way in which nodes affect each other is not only deter- 
mined by their connections, but by logic functions, which 
are generated randomly, simply using lookup tables for each 
node, which take the states of the connecting nodes as in- 
puts, and the state of the node as output. These also remain 
fixed (quenched) during the dynamics of the network. 

We can see that RBNs are a generalization of Boolean cel- 
lular automata (CA) ( |von Neumann, 19661 fWolfram, 1986| 
IWuensche and Lesser, 19 92 1, where the state of each node 
is not affected necessarily by its neighbours, but potentially 
by any node in the network. RBNs with N = K are also 
called random maps. 

The updating of the nodes in classic RBNs is syn- 
chronous: the states of nodes at time f+1 depend on the 
states of nodes at time f, so that all nodes "march in step". 
We will see below that there can be drastic differences if we 



change the updating scheme. 

Usually, an initial random state is chosen for the RBN, 
and the dynamics flow according to the updating functions 
and scheme. Since the state space is finite (2 N ), eventually a 
state will be repeated. Since the dynamics are deterministic, 
this means that the network has reached an attractor. If the 
attractor consists of one state, it is called a point attractor or 
steady state, whereas if it consists of two or more states, it is 
called a cycle attractor or state cycle. The set of states that 
flow towards an attractor is called the attractor basin. An 
example of a RBN with N = 3 and K = 3 is shown in Figure 
□ 

If we try to imagine all possible networks, for each node 
there will be 2 possible functions. And each node has 
N\/(N — K) ! possible ordered combinations for K different 
links. Therefore all the possible networks for given N and K 
will be ( |Harvey and Bossoma ier, 1997| l: 
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Note that many of these will be equivalent, but neverthe- 
less we can see that the space of networks is immense. This 
makes things complicated for statistical studies. There is 
not enough computational power to exhaust all possible net- 
works, so only "representative" properties can be studied. 
We should also mention that there is a very high variance in 
the statistical studies of RBNs. 

However, general properties can be extracted from this 
huge universe of possible networks. 

Order, Chaos, and the Edge 

In RBNs, as well as in many dynamical systems, three 
phases can be distinguished: ordered, chaotic, and critical. 
These phases can be identified with different methods, since 
they have several unique features. 

If we plot the states of a network in a square lattice where 
the state of a node depends topologically on its neighbours, 
and let the dynamics flow, we can easily see which states 
change, and which ones are stable. In other words, we can 
observe how much the network changes. We can colour 
changing states with green, and static ones with red. In the 
ordered phase, we will see that after we select an initial ran- 
dom state, initially many states are changing (green), but 
quicky the dynamics stabilise, and most of the nodes will 
be static (red). There will be only few green "islands", sur- 
rounded by a red "frozen sea". In the chaotic regime, most of 
the states are changing constantly, so we have a green sea of 
changes, typically with red stable islands. The phase tran- 
sition from the ordered to the chaotic regime, also known 
as the "edge of chaos", occurs when the ordered green sea 
breaks into green islands, and the red islands join and perco- 
late through the lattice flKauffman, 2000| pp. 166-167). 



A second feature of these dynamical phases is related to 
"sensitivity to initial conditions", "damage spreading", and 
"robustness to perturbations" which are different ways of 
measuring the stability of a network. We can "mutate", 
"damage" or "perturb" a node of a RBN by flipping its state. 
We can also change a connection between two nodes, or in 
the lookup table of a node. Since nodes affect other nodes, 
we can measure how much a random change affects the rest 
of the network. In other words, we can measure how the 
damage spreads. This can be done by comparing the evolu- 
tion of a "normal" network and a "perturbed" network. In 
the ordered regime, usually the damage does not spread: a 
"perturbed" network "returns" to the same path of the "nor- 
mal" network. This is because changes cannot propagate 
from one green island to another. In the chaotic phase, these 
small changes tend to propagate through the network, mak- 
ing it highly sensitive to perturbations. This is because per- 
turbations can propagate through the percolating green sea 
(Kauffman, 2000, pp. 168-170). This "butterfly effect" is a 
common characteristic of chaotic systems, where small per- 
turbations can cause large consequences, and systems are 
sensitive to their initial conditions. At the edge of chaos, 
changes can propagate, but not necessarily through all the 
network. 

A third feature is the convergence versus divergence of the 
trajectories in state space of the network dynamics. In the or- 
dered phase, similar states tend to converge to the same state. 
In the chaotic regime, similar states tend to diverge. At the 
edge of chaos, nearby states tend to lie on trajectories that 
neither converge nor diverge in state space ( Kauffman, 2000 
p. 171). 

Living systems, or computing systems, need certain sta- 
bility to survive, or to keep information; but also flexibility 
to explore their space of possibilities. This has lead people 
to argue that life and computation occur more naturally at 
the edge of chaos (Langton, 1990), or at the ordered regime 
close to the edge of chaos ( Kauffman, 2000 1. There could be 
ordered or chaotic systems able to perform the same compu- 
tations, but the first ones would need more time, and the sec- 
ond more redundancy to cope with their instabilities. There- 
fore, we can assume that an adequate balance of order and 
chaos is more economic, thus preferable by natural selec- 
tion. 

Phase Transitions in RBNs Very early in the studies of 
RBNs, people realized in simulations that the networks with 
K < 2 were in the ordered regime, and networks with K > 3, 
were in the chaotic regime. In Figure [2] we can appreciate 
characteristic dynamics of RBNs in different phases. 

We can identify phase transitions in RBNs in different 
ways. The main idea is to measure the effect of perturba- 
tions, the sensitivity to initial conditions, or damage spread- 
ing. This is analogous to Lyapunov exponents in continuous 
dynamics. 
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Figure 1: a) Lookup table for the state transitions, b) Wiring diagram: in this case, all nodes affect all nodes, c) State 
space diagram. There is one point attractor (Oil), with one state flowing into it (000), and one cycle attractor of period three 
(111—^110—^101), with three states flowing into it (001, 010, 100). 
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Figure 2: Trajectories through state space of RBNs within 
different phases, N — 32. A square represents the state of a 
node. Initial states at top, time flows downwards, a) ordered, 
K= 1. b) critical, K = 2. c) chaotic, # = 5 



The phase transitions can be statistically or analytically 
obtained. Derrida and Pomeau were the first to determine 
analytically that the critical phase (edge of chaos) was found 
when K = 2 dDerrida and Pomeau, 1986) . They also intro- 
duced two generalizations of the classical model: one where 
they consider nonhomogeneous networks (K is not necessar- 
ily the same for all nodes, so we use as a parameter the mean 
connectivity (K)), and another where the values of lookup 
tables have a probability p of being one (and thus 1 — p of 
being zero). 

The method they used, also known as the Derrida an- 
nealed approximation, takes two random initial configura- 



tions, and measures their overlap. This can be done with 
the normalized Hamming distance (0. Then one time step 
of the dynamics is computed, and the overlap is measured 
again. Then, a new set of rules and connections is cho- 
sen at random. It can be shown that this evolves in a one- 
dimensional map. 
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For p = 0.5, the map converges to the stable point H = 
when K < 2, meaning that different states tend to converge 
(ordered dynamics). At K = 2, this point becomes unsta- 
ble, meaning that for K > 2, different states tend to diverge 
(chaotic dynamics). It can be shown that the curve for criti- 
cal K's depending on the value of p follows the equation (|3ji, 
both for homogeneous and nonhomogeneous normal distri- 
butions of K. The plot of this equation can be seen in Figure 
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Luque and Sole made a simpler analytic determination 
of the phase transitions in networks, where they study 
the damage spreading when single nodes are perturbed 
( |Luque and Sole, 199 7b). They consider trees of nodes that 
can affect the state of other nodes in time. As a node has 
more connections, there will be an increase in the proba- 
bility that a damage in a single node (0— >1 or 1— >0) will 
percolate through the network. 
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Figure 3: Phase diagram for the classical model, reprinted 
from ( |Aldana, 20 03 1 with permission from Elsevier 



Let us focus only in one node i at time t, and a node j 
of the several i can affect at time t + 1. There is a prob- 
ability p that j will be one, and a damage in i will mod- 
ify j towards one with probability 1 — p. The complemen- 
tary case is the same. Now, for K nodes, we could expect 
that at least one change will occur if (K)2p{ \ — p) > 1, 
which leads to Q. This method can be also used for other 
types of networks. Luque and Sole later used the concept of 
Boolean derivative to define Lyapunov exponents in RBNs 
( |Luque and Sole, 2000) : 



X = \og[2p{\-p)K] 



(4) 



where A, < represents the ordered phase, A, > the 
chaotic phase, and A, = the critical phase. 

Statistical studies confirm these analytical results. How- 
ever, it seems that, in practice, the size of the network 
can play a role in the phase transitions (Gershenson, 2004a ; 
IGershenson, 200 4b I. We have seen that for large RBNs the 
phase transition (measured with the average of differences of 
normalized Hamming distances at t — > °° and t = 0, for min- 
imum initial distances of 1 /N) is given for K shifted towards 
one, and for very small networks for K shifted towards three. 
This is probably because for large networks, there is a higher 
probability that a subnetwork will be generating more noise 
than "average", thus propagating damage. 

In practice, there are very high standard deviations in 
RBN studies, which are normal in this type of systems 
(Mitc hell et al., 19931 - We can clearly find (or design) or- 
dered networks for K 3> 3 but on average, statistical studies 
confirm the analytical ones. Therefore, when we generate a 
random network, there are high probabilities that it will be 
in a specific regime according to K and p. 



Explorations of the Classical Model 

There have been several explorations of differ- 
ent properties of RBNs, e.g. ( fWuensche, 1997| 
Alda na-Gonzalez et al., 2003V One can measure, for 
example, the number and length of attractors, the sizes 
and distributions of their basins, and how these depend 
on different parameters of RBNs, such as N, K, p, or the 
topology. 

The structure of the nodes is very important for the dy- 
namics of RBNs. The descendants of a node are the nodes 
that it affects, while the ancestors of a node are those that 
affect it. To have cycle attractors, i.e. of period greater 
than one, there should be at least one node that will be its 
own ancestor. A circuit of auto-activating nodes is called a 
linkage loop, and when there is no feedback, linkage trees 
are formed. Note that loops spread activation through trees, 
but not vice versa. The relevant elements of a network are 
those nodes that form linkage loops, and do not have con- 
stant functions, for these cause instabilities in the network, 
which might or not propagate. Note that as there are more 
connections in a network (higher K), the probability of hav- 
ing loops increases. Therefore, finding less stable dynamics 
for high values of K is natural. 

We should not confuse the node diagrams with the state 
space diagrams. These show the dynamic trajectories of 
states of the network, while the first show the relations of the 
network elements. Classic RBNs are dissipative systems: a 
state can have only one successor, since the dynamics are 
deterministic, but more than one predecessor or pre-image 
can flow into a single state. The in-degree of a state is the 
number of predecessors it has. States without predecessors 
are called garden- of -Eden states. In this way, the dynamics 
flow from garden-of-Eden states, converging towards attrac- 
tors. The time it takes to reach an attractor is called transient 
time. 

Attractor Lengths There have been analytic solutions of 
RBN s foxK = 1 (|Flyvbjerg and Kjaer, 19 88}, and for K = N 
( |Derrida and Fly vbjerg, 1987 1. The challenging problem of 
finding a general analytic solution is still open. Some sta- 
tistical studies have matched the analytic solutions for the 
special cases. People have observed the following, consid- 
ering p = 0.5 (|Kauffman, 1993| IBastolla and Parisi, 1998] 
Aldana-Go nzalez et al., 2003} : 

For K = 1, the probability of having long attractors de- 
creases exponentially, and the average number of cycles 
seems to be independent of N ( Bastolla and Parisi, 1998 1. 
The median lengths of state cycles are of order ^/N /2. 

For K >N, the average length of attractors and the tran- 
sient times required to reach them grow exponentially. This 
restricts numerical investigations to small networks. The 
typical cycle length grows proportional to 2 N I 2 . 

For K = 2, at the critical phase, both the typical at- 
tractor lengths and the average number of attractors grow 



algebraically with N. However, the precise dependence 
of N is a matter of dispute. People long believed that 
the average number of attractors and their length was 
proportional to VN, (|kauffman ,T969| |Kauffman, 19931 
IBastolla and Parisi, 1998| . This result was very attractive, 
because the number of cell types and cell replication times 
for different organisms seem to scale also as the square 
root of genes for different species, although the precise 
number of genes of organisms keeps on changing. How- 
ever, Bilke and Sjunnesson did a full exploration of net- 
works, "decimating" irrelevant variables, and found that 
there is a linear dependence of number of attractors depend- 
ing on N (Bilke and Sjunnesson, 2002 1. This linear depen- 
dence has been confirmed in other complete statistical stud- 
ies ( Gersh enson, 200 2 Gersh enson et al., 20031 . The differ- 
ence seems to lie on the bias caused by undersampling the 
state space. 

Since there are 2 N states, full statistical studies are possi- 
ble only for very small networks. For example, for N = 20 
there are more than a billion initial states. Therefore, these 
studies either concentrate on small networks, or take into 
account only very few initial states. In the first case, some 
properties of large networks will not be observed, whereas in 
the second, some attractors would not be found, especially 
if their basins consist of very few states. 

More research is needed in this direction. One argument 
could be that even when there would be potentially more at- 
tractors than the ones found by limited sampling, in practice 
one would obtain the same result, since nature does not ex- 
haust all possible configurations (e.g. there could be more 
cell types for a given number of genes, but developing into 
them is impossible, i.e. their basins are very small). Also, if 
we could get the general analytical solution, it could be that 
it would not match statistical studies, since a bias should be 
expected in the statistical sampling, due to very high stan- 
dard deviations. However, for different purposes we might 
be more interested in the practical than the theoretical re- 
sults, or vice versa. 

Convergence One can measure the convergence of states 
with different parameters. One of them is the G-density, 
which is the density of garden-of-Eden states. Another is 
the in-degree frequency distribution, which can be plotted 
as a histogram (Wuensche, 1998). These measures reveal 
features at different phases. 

At the ordered phase, there is a very high G-density, and 
high in-degree frequency. This leads to a high convergence, 
and very short transient times. The basins of attraction are 
very compact, with many states flowing into few states. 

At the critical phase, the in-degree distribution approx- 
imates a power-law, i.e. there are few states with high 
in-degree, and many states with low in-degree. There is 
medium convergence. 

At the chaotic phase, there are a relatively lower G- 



density, and a high frequency of low in-degrees. The basins 
of attraction are very elongated, with few states flowing into 
other states. This makes average transient times very long, 
and in some cases infinite in practice. Therefore, there is low 
convergence. 

Other parameters that can be useful for measuring 
convergence include Walker's "internal homogene- 
ity" (Walker and Ashby, 1966 1, Langton's A. param- 
eter (|LangtonrT990X and Wuensche's Z parameter 
dWuensche, 1999> . The latter one, together with the "input- 
entropy variance", can be also used to automatically classify 
rules of CA into ordered, complex (critical), and chaotic 
( |Wuensche, 1999) . This is useful, since "interesting" 
behaviour in CA tends to occur within complex rule space. 

Multi- Valued Networks 

There have been some extensions to the Boolean idealization 
of classical RBNs, namely where nodes can take more than 
two values. 

Sole, Luque, and Kauffman, and more recently Luque 
and Ballesteros have studied such multi-valued networks, 
and calculated their phase transitions ( |Sole et al., 2000| 
|Luque and Ballesteros, 2 004 1. For the special case where 
only two states are allowed, the results of Derrida are re- 
covered. 

In nature, the components of certain systems exhibit a be- 
haviour that is better described with more than two states. 
Particular models should go beyond the boolean idealiza- 
tion. However, for theoretical purposes, we could combine 
several boolean nodes to act as a multi-valued one (codify- 
ing in base two its state). 

Topologies 

Many systems have been found to have a scale-free topol- 
ogy. It seems to be a persistent feature of complex networks 
(Barabasi, 2002): The Internet, molecular and genetic net- 
works, social networks, technology graphs, language net- 
works, food webs... they all share similar topological fea- 
tures: they have few elements with many links, and many 
elements with few links. This distribution seems to have 
several adaptive advantages, thought it is still not very well 
understood. 

However, most RBNs which have been studied have ho- 
mogeneous or normal topologies. Oosawa and Savageau 
studied the effects of topology in the properties of RBNs 
( |Oosawa and Savageau, 2002 1. Their results showed that the 
topology can change drastically these properties. Networks 
with the more uniform rank distributions exhibit more and 
longer attractors and less entropy and mutual information 
(less correlation in their expression patterns), whereas more 
skewed topologies exhibit less and shorter attractors and 
more entropy and mutual information. A topology based on 
E. coli, which is scale-free, balances the parameters to avoid 
the disadvantages of the extreme topologies. 
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Figure 4: Phase diagram for the scale-free model, reprinted 
from flAldana, 20 03 1 with permission from Elsevier 

Having this in mind, Aldana studied many properties of 
RBNs with scale-free topology ( Aldana, 2003 1. The connec- 
tivity of scale-free RBNs can be generated randomly using 
the probability distribution P(k) = (y) ~ 1 , where y > 1 
and i^(y) = J2T=i^ ^ * s trie Ri emann Zeta function. In this 
way, every node has at least one connection, but there are 
few ones with many connections. The properties of the net- 
work are no longer determined by the average connectivity, 
but by the exponent y. 

Following Derrida's method, Aldana found that the criti- 
cal value of the exponent y c where the phase transition from 
order to chaos occurs is determined by the transcendental 
equation: 

2P(1-P)^-^ = 1 (5) 

The values for which (JSJl is satisfied are plotted in Figure 
|4] We can see that y c e [2,2.5] for any value of p. The 
maximum value of y 6 . w 2.47875 is reached when p = 0.5. 

The network properties at each phase (e.g. number and 
length of attractors, transient times) are analogous to the 
ones obtained with homogeneous RBNs. 

An important result is that evolvability has more space in 
scale-free networks, since these can adapt even in the or- 
dered regime, where changes in well-connected elements do 
propagate through the network. However, experimental ev- 
idence shows that most biological networks are scale free 
with exponent 2 < y < 2.5 ( Aldan a72003t , i.e. "at the edge 
of chaos". The advantages of scale-free topologies are be- 
ginning to become evident, although they have been not em- 
braced by most researchers of RBNs. 

RBN Control 

RBNs usually do not consider external inputs. However, real 
systems such as genetic networks can be influenced by ex- 



ternal signals, such as molecular clocks related to sunlight. 

Methods of chaos control have been successfully 
applied to chaotic RBNs (Luque and Sole, 1997a 
Luq ue and Sole, 1998| |BaUestero^andLTqueT 2002 ). 
The main idea is to use a periodic function to drive a very 
chaotic network into a stable pattern. If a periodic function 
determines the states of some nodes at some time, these 
will have a regularity that can spread through the rest of the 
network, developing into a global periodic pattern. A high 
percentage of nodes should be controlled to achieve this. 
However, once we control a small chaotic network, we can 
use this to control a larger chaotic network, and this one 
to control an even larger one, and so on. This shows that 
it is possible to design chaotic networks controlled by few 
external signals to force them into regular behaviour. 

Different Updating Schemes 

Kauffman has argued that the small average number of at- 
tractors found in RBNs compared with the number of possi- 
ble states can account for the number of cell types and cell 
replication time in organisms compared with their number 
of genes ( |Kauffman, 1969||Kauffman, 1993) . In those days, 
about eighty thousand genes were thought to conform the 
human genome. Therefore, if the genome is seen as a RBN 
close to the edge of chaos, the expected number of attractors 
would be less than three hundred, matching the observed 
number of cell types in humans. However, there are many 
drawbacks to this calculation. The Boolean idealization has 
been roughly accepted, since multi-valued networks have 
shown similar results. Now that the human genome has been 
mostly sequenced, it seems to consist of less than thirty five 
thousand genes. The topology seems to be scale-free. There 
is certain amount of junk or structural DNA, without func- 
tionality. Many functions seem to be biassed (p ^ 0.5). But 
the heaviest argument has been the following: genes do not 
march in step. Genes do not change their states all at the 
same moment, but some do it earlier than others. There was 
no argument for the synchronicity in RBNs. In the next sec- 
tions we will review research made related to this criticism. 

Asynchronous RBNs 

Harvey and Bossomaier introduced the criticism to the syn- 
chronicity of classic RBNs ( H arvey and Bossomaier, 1997| l. 
It was well known that asynchronicity could change dras- 
tically the dynamics of a synchronous system, such as the 
prisoner's dilemma ( |Huberman and Glance, 1993} or Con- 
way's game of life (Bersini and Detours, 1994), and they did 
a similar thing for RBNs. Instead of updating the nodes 
synchronously, they defined asynchronous RBNs (ARBNs), 
where a node is picked up at random, and updated. We 
have to notice that ARBNs are not only asynchronous, but 
also non-deterministic. This destroys the cycle attractors of 
classical RBNs (CRBNs), since it is very difficult that a se- 
quence of states will be repeated with a non-deterministic 
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Figure 5: State space diagram of an ARBN, using the lookup 
table and wiring of Figure^ Dotted arcs indicate stochastic 
transitions. There is one point attractor (011) and one loose 
attractor (100,101,110,111). The state (0,0,0) could lead to 
either attractor, via (010) or (001), respectively; i.e. it is in 
the basin of both attractors 



updating. Point attractors still appear in ARBNs. There are 
also "loose" attractors, which can be seen as a subset of the 
state space which "traps" the dynamics after some time (all 
the possible states would rarely be visited). 

The behaviour of ARBNs changes drastically from the 
one presented by CRBNs. Not only cycle-attractors disap- 
pear, but their basins can change. Also, states in ARBNs can 
be in more than one basin of attraction: they have the po- 
tentiality to fall into different attractors depending on which 
nodes are updated. An example of an ARBN can be appre- 
ciated in Figure |5] 

Harvey and Bossomaier found that in theory, there is on 
average exactly one point attractor for any ARBN family. 
However, since we do not exhaust all possible networks, a 
skewed distribution of the point attractors is made evident: 
for some values ARBNs have more skewed distributions, i.e. 
few networks with many point attractors, several with none, 
so the averages tend to be lesser than one (except for the 
special case K = 0, where it is exactly one), reaching a min- 
imum forK = 3 (Gershenson, 2002 1. It seems that there are 
no loose attractors for K = 1, and the probability of having 
one increases with K. 

The properties of ARBNs, being so different to the ones 
of RBNs, casted a doubt on the validity of CRBNs as models 
of genetic regulatory networks. 

Rhythmic Asynchronous RBNs If ARBNs were to 
model biological systems, how could they have any rhythm? 
An artificial way of solving this could be to imple- 
ment a CRBN in an ARBN using Nehaniv's method 
(Nehaniv, 2002 1, but such a network would be very unre- 



alistic. Another option could be to introduce proportionally 
large time delays in each node, so that each node could be 
updated only when most probably all the other nodes would 
have been updated (Klemm and Bornholdt, 2003 1. How- 
ever, this is in a way a disguised synchronicity, and not more 
realistic than it. 

Trying to find a solution, Di Paolo used genetic algorithms 
to explore the space of possible ARBNs, and found networks 
that do have rhythmic behaviour ( Di Paolo, 2 00 1 1 . Recently, 
Rholfshagen and Di Paolo analysed the topology of such 
rhythmic ARBNs (Rholfshagen a ndDi Paolo, 2004) . They 
found that invariably this consists of a "ring" of nodes that 
affect only one to the next, i.e. a linkage loop. The rest of 
the nodes affect nodes only outside the ring, so that activa- 
tion can spread only outwards of the ring. In other words, 
there is a single linkage loop in the network, and the dynam- 
ics of this propagate only towards linkage trees. The num- 
ber of nodes in the ring determines the average period of the 
rhythm in epochs (time steps x AO- This is because only one 
node can change the guiding dynamics of the network at a 
time, and any node is updated on average once an epoch. 
These results are very promising, but there are open tasks, 
such as the exploration of rhythmic ARBNs with more than 
one rhythmic attractor. 

Deterministic Asynchronous RBNs 

I agree with the criticisms to the synchronous assumption: 
genes do not march in step. But I do not believe that they 
are random. Having this in mind, I proposed Deterministic 
Asynchronous RBNs (DARBNs) (Gers henson, 2002| l. 

For having nodes that do not update simultaneously, we 
can introduce two parameters per node, P, and Qj (Pj,Qi € 
N,Pi > Qi). All P('s and Qi's are generated randomly with 
maxima P max and Q max , and remain fixed. A node i will be 
updated when the modulus of the time step over Pi equals 
Qi. Like this, P, can be seen as the period of the node up- 
date, i.e. the number of time steps that will pass between 
updated, and Q, can be seen as the translation of the node 
update. If at a certain time step more than one node fulfills 
its updating condition, then the nodes will be updated in an 
arbitrary order (e.g. from left to right), and the whole net- 
work will be updated with each node. 

We can also define Deterministic Generalized Asyn- 
chronous RBNs (DGARBNs), where if more than one 
node fulfills its updating condition, all of these nodes will 
be updated synchronously. This makes DGARBNs semi- 
synchronous. 

For completeness, we also defined Generalized Asyn- 
chronous RBNs (GARBNs), which are like ARBNs, only 
that at each time step, some nodes are randomly selected, 
and these are updated synchronously. In this way, GARBNs 
are semi-synchronous, but non-deterministic. We can see 
examples of these RBNs in Figure[6] 

We found out that all types of networks have the same 
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Figure 6: State space diagrams of a DARBN, a DGARBN, and a GARBN, using the lookup table and wiring of Figure 
H For the deterministic networks, P — {1,1,2} and Q — {0,0,0}. Numbers near arcs indicate the transitions which will 
occur at time modulus two. There is one point attractor (011) and one cycle attractor, but different than the one of a CRBN 
(100i — > lllo — * llli — > 110o), but the basins are different. For the GARBN, there is the same point attractor and loose 
attractor than for the ARBN, but there are more possible stochastic transitions, i.e. changes can potentially propagate faster 



point attractors (Gershenson, 2002 1. This is because no mat- 
ter which node is updated, the state will not change. There- 
fore, the updating scheme does not affect point attractors. 
However, the attractor basins do change, and other types of 
attractor. DARBNs and DGARBNs have properties much 
closer to the ones of CRBNs ( Gershenson, 2002 1. Since they 
are deterministic, cycle attractors are present. The number 
of attractors, like in CRBNs, increases linearly with N, but 
more slowly, meaning that there are even fewer attractors for 
a high possible number of states 1 . Also, the percentage of 
states in attractors is reduced exponentially, as with CRBNs, 
but even faster. These results imply that DGARBNs and 
DARBNs can perform even more complexity reduction than 
CRBNs. We also concluded that the difference of CRBNs 
and ARBNs lies more in non-determinism than in asyn- 
chronicity. 

Moreover, we proposed a method for mapping any deter- 
ministic asynchronous RBN into a CRBN. The main idea is 
to introduce new nodes, connected to every other previous 
node, which codify in base two the maximum period of the 
DARBN. This is the least common multiple of all Pi's. 

Asynchronicity and Feedback We have to mention that 
Thomas developed much earlier an asynchronous model 
of RBNs using delays, which could be both determinis- 
tic or stochastic, depending on the certainty of the delays 
( |Thomas, 19731 IThomas, 19781 IThomas, 1991) . However, 
these models were proposed and used mainly for the analysis 
of precise networks, their circuits, and feedback loops, and 
to my knowledge they have been not used for analytical or 
statistical studies of ensembles ("families") of networks. An 



1 ARBNs and GARBNs have also a linear increase in the aver- 
age numbe r of attractors, wh en we consider loose attractors in the 
statistics l Gershenson, 2004b I. These RBNs have less attractors 
than deterministic RBNs, but their sizes are much larger. 



interesting finding of Thomas and coworkers was the follow- 
ing: a positive feedback loop (direct or indirect autocataly- 
sis) in a network implies the choice of two stable states, i.e. 
point attractors. This gives the property of multistationarity . 
On the other hand, a negative feedback loop implies peri- 
odic behaviour, i.e. point or cycle attractors. This can be 
described as homeostasis. The combinations of positive and 
negative feedback loops give RBNs a plethora of possible 
behaviours, many of which are found in living systems. 

Mixed-context RBNs 

We can see the sets P and Q, consisting of all P,'s and Qi's 
respectively, as the context of a network. This is because 
external factors, such as temperature, can change the updat- 
ing periods of elements of a system. Like this, the same 
DGARBN can have different behaviours in different con- 
texts, i.e. for different P and Q. 

Having this in mind, we can introduce non- 
determinism in contextual RBNs in a very specific 
way (Gershen son et al., 2003| >. For a given DGARBN, 
we can have M "pure" contexts. Then, each R time steps 
we select randomly one of these contexts, and use it in 
the network. Thus, we have defined Mixed-context RBNs 
(MxRBNs). An example is shown in Figure[7] 

We should note that the non-determinism in MxRBNs is 
introduced in a very controlled fashion. In GARBNs we 
have N "coin flips" per time step (selecting which nodes 
will be updated). In ARBNs we have one coin flip per time 
step (selecting which node will be updated). In MxRBNs 
we have one coin flip per R time steps (selecting which con- 
text will be used). The higher the value of R and the lower 
number of M contexts, the less stochasticity there will be. 

MxRBNs have very similar number of attractors than 
ARBNs and GARBNs (considering loose attractors), and 
also their number increases linearly and slowly with N. 
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Figure 7: State space diagram of a MxRBN, using the 
lookup table and wiring of Figure [2 M = 2 contexts. 
Pi = {1, 1,2},? 2 = {2, 1,1} and Q Y = Q 2 = {0,0,0}. There- 
fore, the dynamics are deterministic and context indepen- 
dent when the modulus of time over two is zero, and depend 
on the context when it is one. There is one point attractor 
(0 1 1 ) and one loose attractor ( 1 00, 1 1 , 1 1 0, 1 1 1 ) 

However, they can perform significantly more complexity re- 
duction ( Gersh enson, 200 4b ), even more than CRBNs, but 
not as much as DGARBNs or DARBNs. This means that 
very few states, from all their possible states, lie in their at- 
tractors. 

We can see that the context in RBNs allows a high in- 
crease in complexity reduction, since it allows information 
to be "thrown" into the context. We can restate this with 
the following argument. Any updating scheme can perform 
in theory any computation. The more non-determinism we 
introduce, the more redundancy in the network we need to 
perform the computation. We require more elements. On 
the other hand, contextual RBNs can exploit information in 
their contexts, which CRBNs need to codify explicitly. In 
terms of number of elements required for a computation, 
less will be required in a DGARBN. Then the ascending or- 
der would be of DARBNs, MxRBNs, CRBNs, ARBNs, and 
finally G ARBNs. 

Classification of RBNs 

We can classify different types of RBNs according to their 
updating scheme (Gers henson, 2002) . All RBNs are dis- 
crete dynamical networks (DDNs) (Wuensche, 1997 1, since 
they have discrete time, states, and values. The most gen- 
eral type of RBNs are GARBNs, since all of the others 
can be seen as special cases of them. If, on one hand, we 
make them deterministic with a context conformed by P 
and Q, then we would have DGARBNs. MxRBNs would 
be more general than DGARBNs, since they have M con- 
texts. DGARBNs are a special case, when M = 1 or R — > 



DDK 




Figure 8: Classification of random Boolean networks, ac- 
cording to their updating scheme 

oo. Random maps would be a special case, of DGARBNs 
when P = 1, Q = 0, and N = K. Random maps can 
simulate with redundancy any CRBN, but not vice versa, 
so the latter can be seen as a special case of the former. 
Boolean CA are special cases of CRBNs, where the con- 
nectivity and functionality are the same, i.e. symmetri- 
cal, for all nodes. From GARBNs, on the other hand, we 
can limit to the update of only one node at a time, and 
we will have ARBNs ( |Harvey and Bossomaier, 1997) . If 
we make them deterministic with a context, then we have 
DARBNs. These and DGARBNs overlap on the special 
cases when one and only one node is updated at a time 
step. There are special ARBNs with rhythmic and non- 
rhythmic attractors ( Di Paolo, 2001| l, and probably there are 
also GARBNs with these types of attractors. We should 
note that particular instantiations of Thomas' asynchronous 
RBNs ( IThomas, 1973l|Thomas, 1978||Thomas, 1991| l could 
be seen as ARBNs or DARBNs, depending on the certainty 
of the delays. We can see a diagram of the proposed classi- 
fication in Figure|8] 

We could add a third dimension to Figure[8] indicating the 
number of allowed states per node, to include multi-valued 
networks. Special cases of the different RBNs presented 
can be also interesting, such as networks with scale-free or 
another particular topology, or CA with different updating 
schemes. We can see that DDNs offer a rich variety of mod- 
els, many of which have not been yet studied. 

The more general a type of RBN is, the more trajectories 
in state space it can have. Note that since we only change the 
updating scheme, we can easily compare the behaviour of 
the same network, i.e. rules and connectivity, with different 
schemes. We have seen that the dynamics and the basins 
of attraction can change considerably as we switch updating 
scheme, although in extreme cases, such as K — or when 
the number of attractors equals N, different network types do 
have the same behaviour. An example can be seen in Figure 
[9] However, we have recently found that, even when the 
precise stability of RBNs depends on the updating scheme, 
the phase transition seems to be very similar independently 
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Figure 9: Dynamics of the same RBN under different updating schemes, with N = 32, K = 2,p = 0.5, for the same initial state. 
For contextual RBNs, maxP = 5,maxQ = 4. For MxRBN, M = 2,R = 10. Few states are frozen (e.g. second from left to right 
is always black, the fourth is always white), but dynamics change drastically with updating scheme. Note that in spite of being 
non-deterministic, the dynamics of the MxRBN are visually more similar to the ones of DARBN or DGARBN than ARBN or 
GARBN. 



of the updating scheme (Gershens on, 2004at . 

Applications 

The generality of RBNs makes them applicable in a wide 
range of domains. In this section we review some of the 
main areas where RBNs have been applied. 

Genetic Regulatory Networks 

Most cells in multicellular organisms have the same ge- 
netic information, although the genes that are expressed or 
"on" change at every moment. Genes interact with each 
other via proteins, but there are so many of them, that ge- 
netic regulation is not fully understood. RBNs were orig- 
inally proposed to tackle this problem (Kauff man, 1969) . 
They have been used not only to explore their generic prop- 
erties, but also to analyse and predict genomic interac- 
tions (|Somogyi and Sni eg oski, 1996| ISomogyi et al.,"T997| 
|D'haeseleer et al., 1998| l. Since the genomic data is incom- 
plete, and certain knowledge about real genetic networks 
is required for disease treatment, probabilistic boolean net- 
works (PBNs) were introduced by Shmulevich and cowork- 
ers. These are useful for inferring possible gene functional- 
ity from incomplete data ( Sh mulevich et al., 2 002 1. 

There has been experimental evidence for Kauff- 
man's interpretation of cell types as attractors of RBNs 



( |Huang and Ingber, 2000| l, which encourages more RBN re- 
search. It seems that many genes do not play a role in the 
type of a cell. They certainly might have other functions, 
for example related to the cell's metabolism. Nevertheless, 
there is a very strong correlation for some genes as a cell 
type is mechanically forced, which shows that the activation 
patterns of a subset of genes can indicate the type of a cell. 

There have been also other models of genetic reg- 
ulatory networks which include continuos states 
( |Glass and Kauffman, 19731 |Kappler et al., 2002| i. Us- 
ing differential equations in which gene interactions are 
incorporated as logical functions, there is no need for a 
clock to calculate the dynamics. 

Evolution and Computation 

RBNs are very tempting models for studying evolution, 
since we do not assume any functionality. It is now known 
that evolvability is expected at the edge of chaos, where 
small changes do not destroy previous functionality, but can 
explore their space of possibilities incrementally. RBNs also 
present naturally modules, a desired feature in evolvable sys- 
tems. 

Fernandez and Sole have studied issues related to the 
evolvability of networks, such as robustness, redundancy, 
degeneracy, and modularity (Ferna ndez and Sole, 2 004). 
This is a very abstract study of the requirements of life, 



since one way we can distinguish physical from biolog- 
ical systems is to note that the latter perform computa- 
tions (Hopfield , 1994) . By understanding how a network 
can evolve to perform certain computations, we are answer- 
ing questions on how complex organisms evolved on Earth. 
There have been some studies of evolution of RBNs us- 
ing genetic algorithms with promising results ( |Stern, 19991 
Lemk eetal., 2001) . 

Related work is the one carried out in the area of evolvable 
hardware, where complex logical circuits are evolved in re- 
configurable hardware (Thompson, 199SJ. Also, research in 
RBNs could provide valuable feedback for evolving logical 
circuits. 

Other areas 

RBNs have been used in several other areas, such as neural 
networks ( |Huepe- Minoletti and Aldana- Gonzalez, 2002) , 
social modelling j jSheTHng, 1971), robo tics 
( |Quick et al., 2003) , and music generation ( |Dorin, 200"0Y 

Finally, RBNs are interesting mathematical objects by 
themselves. Since cellular automata are special cases of 
RBNs, we can find many more applications there, e.g. in 
percolation theory (Stauff er, 1985| l. 

Tools 

There are several software applications available for the ex- 
ploration of different properties of RBNs: 

DDLab. Developed by Andy Wuensche, Discrete Dy- 
namics Lab is probably the most powerful tool for studying 
discrete dynamical networks: synchronous RBNs and CA, 
including multi-valued networks. One can observe the dy- 
namics of the networks, explore their basins of attraction. 
It includes a wide variety of measures, data, analysis and 
statistics. It is very well documented, and runs on most plat- 
forms. It is available at http://www.ddlab.com 

RBN Toolbox for Matlab was developed by Christian 
Schwarzer and Christof Teuscher. It is useful for simulat- 
ing and visualizing RBNs. It includes different updating 
schemes, statistical functions, and other features. It is avail- 
able at http://www.teuscher.ch/rbntoolbox 

RBNLab was developed by the author. It can graphi- 
cally simulate the dynamics of RBNs with different updating 
schemes. It can find attractors (point, cycle, and loose), and 
generate different statistics. It is implemented in Java, so it 
runs on most platforms. The source code and the program 
are available at http://rbn.sourceforge.net 

BN/PBN Toolbox for Matlab is maintained by Harri 
Lahdesmaki and Ilya Shmulevich. It can be used to work 
with CRBNs and Probabilistic Boolean Networks. It in- 
cludes functions for simulating the network dynamics, com- 
puting network statistics (numbers and sizes of attractors, 
basins, transient lengths, Derrida curves, percolation on 2-D 
lattices, influence matrices), computing state transition ma- 
trices and obtaining stationary distributions, inferring net- 



works from data, generating random networks and func- 
tions, visualization and printing, intervention, and mem- 
bership testing of Boolean functions. It is available at 
http://www2.mdanderson.org/app/ilya/PBN/PBN.htm 

Future Lines of Research 

There are many promising lines to be followed in RBN 
research. The ensemble approach is promising for un- 
derstanding and predicting properties of cells and organ- 
isms (Kauffman, 2004 1. The use of RBNs for data min- 
ing and genetic network analysis is being propagated 
(D'haesel eeretal., 19981 IShmulevich et al., 2002) . Be- 
cause of their generality, RBNs are also interesting for 
studying evolvability and adaptability at an abstract level. 
Generalizations, combinations, and refinements of the dif- 
ferent types of RBNs presented are also worth pursuing, e.g. 
ensemble studies on Thomas' ARBNs or PBNs, scale-free 
multi-valued DGARBNs, etc. Also, the challenging prob- 
lem of finding analytical solutions for CRBNs is still evad- 
ing us. In short, there is plenty of RBN research to do in the 
years to come. 

Conclusions 

This tutorial was a brief introduction to some of the research 
related to random Boolean networks in the past decades. 
There has been much more work that could not be included 
because of different constrains. However, the readers are 
invited to deepen their knowledge in RBNs following the 
references of this tutorial. 

There will be many arguments on which will be the "best" 
model for different purposes and phenomena. CRBNs may 
be not very close to reality, but full ARBNs are even far- 
ther, since genes are not updated in a fully stochastic man- 
ner 2 . Models such as DGARBNs or MxRBNs seem to be 
more realistic, although they are more complicated. For 
ensemble studies, for simplicity, CRBNs can be justifiable 
( |Gershenson, 20 04b I. For modelling particular genetic net- 
works, other models might be more realistic and useful, such 
as PBNs, or the models of Thomas, Glass, or Somogyi. 
However, their complexity would make difficult, but still in- 
teresting, ensemble studies on them. We can say that differ- 
ent models will be more suitable for different purposes, and 
most of them are worth exploring. 

We can say that RBNs are very inviting, because of their 
generality: one can model at a very abstract level many phe- 
nomena, and study generic properties of networks indepen- 
dently of their functionality. They have also raised important 

2 We have found out that ARBNs and GARBNs perform 
less complexity redu ction than deterministic RBNs or MxRBNs 
I Gershenson, 2004a l. This indicates that determinism or quasi- 
determinism is a desirable property of natural systems. How could 
this evolve, is a diffe rent question, a nd studies such as the ones 
on rhythmic ARBNs I Di Pa olo, 200 1) are providing interesting an- 
swers. 



questions related to the differences between theory and prac- 
tice, since some analytical or statistical results do not match 
each other. This can be explained with skewed distributions 
and very high standard deviations found in RBN statistics. 
However, this brings deeper philosophical questions: how 
much should we care about theory, and how much should 
we care about practice? Which one will help us understand 
better the phenomena we try to study? It seems that a care- 
ful balance of both is required. However, it is unknown how 
this balance might be like. 
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